#OS level tools
import os
import time
import itertools
from collections import defaultdict
from glob import glob
import psutil
from functools import partial
import re
from __future__ import print_function
#array and data structure
import numpy as np
import pandas as pd
#import seaborn as sb
#Ipython display and widgets
#import ipywidgets as widgets
from IPython.display import Image, HTML, display
from IPython.display import Markdown as md
#from ipywidgets import interact_manual
#holoviews and plotting
import holoviews as hv
import datashader as ds
from holoviews.operation.datashader import aggregate, shade, datashade
from bokeh.models import HoverTool
#from holoviews.operation import decimate
#dask parallelization
import dask.dataframe as dd
from dask import compute, delayed
import dask.threaded
import dask.multiprocessing
#tsne
from MulticoreTSNE import MulticoreTSNE as TSNE
tsne = TSNE(n_jobs=24)
#color assignment
cmap_all=['white','white']
cmap_parent=['black','grey']
cmap_pop=(['darkgreen','lightgreen'], ['darkorange','yellow'], ['purple','blueviolet'], ['darkblue','lightblue'], ['indianred','red'])
for i in range(5):
cmap_pop=cmap_pop+cmap_pop
background = '#D3D3D3'
#export path assignment
#scratch_path='/scratch/'+os.environ['USER']+'/'+os.environ['SLURM_JOBID']
scratch_path="."
export_path="PNG"
png_path="PNG"
try:
os.makedirs(export_path)
except OSError as e:
if e.errno != os.errno.EEXIST:
raise
#export = partial(export_image, export_path=export_path, background=background)
hv.notebook_extension('bokeh')
#display(HTML("<style>.container { width:100% !important; overflow-x: auto;white-space: nowrap;}</style>"))
hv.opts("RGB [toolbar=None, width=400, height=400, bgcolor='#D3D3D3', fontsize={'title':8, 'xlabel':8, 'ylabel':8, 'ticks':3}]")
import ipywidgets as widgets
box_layout = widgets.Layout(overflow_x='scroll',
border='3px solid black',
height='',
flex_flow='column',
display='flex')
row_layout = widgets.Layout(min_width='32000px')
hv.opts("RGB [toolbar=None, width=400, height=400, bgcolor='#D3D3D3', fontsize={'title':15, 'xlabel':10, 'ylabel':10, 'ticks':5}]")
def config_objects(s):
try:
with open(s) as config_file:
config_file.seek(0)
gates={}
for line in config_file:
phenoType=""
line = line.strip()
gate = line.split("\t")
if len(gate)==12:
phenoType=gate[11]
gates.update({"pop"+str(gate[0]):[int(gate[0]), int(gate[1]), int(gate[2]), int(gate[3]), int(gate[4]), int(gate[5]), int(gate[6]), int(gate[7]), int(gate[8]), int(gate[9]), int(gate[10]), phenoType]})
return gates
except:
raise Exception("Error parsing configuration file")
def config_summary(s, h):
try:
with open(s) as config_file:
config_file.seek(0)
gates={}
for line in config_file:
phenoType=""
line = line.strip()
gate = line.split("\t")
xmarker=str(h[int(gate[1])-1])
ymarker=str(h[int(gate[2])-1])
startx=int((float(gate[3])/200)*4096)
starty=int((float(gate[5])/200)*4096)
endx=int((float(gate[4])/200)*4096)
endy=int((float(gate[6])/200)*4096)
parent="pop"+gate[7]
if len(gate)==12:
phenoType=gate[11]
gates.update({"pop"+str(gate[0]):[int(gate[0]), parent, xmarker, ymarker, phenoType, startx, endx, starty, endy]})
return gates
except:
raise Exception("Error parsing configuration file")
_nsre = re.compile('([0-9]+)')
def natural_sort_key(s):
return [int(text) if text.isdigit() else text.lower()
for text in re.split(_nsre, s)]
def natural_sort(l):
#https://stackoverflow.com/a/4836734/846892
convert = lambda text: int(text) if text.isdigit() else text.lower()
alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ]
return sorted(l, key = alphanum_key)
def label_color (pops, row):
eventcolor=0
for i, pop in enumerate(pops):
if row[pop]==0:
eventcolor=i+1
return eventcolor
def label_color2 (pops, row):
eventcolor="base"
for i, pop in enumerate(pops):
if row[pop]==0:
eventcolor=pop
return eventcolor
def parseCohort(s):
cohort_file=open(s)
return
def parseDataFrame(s):
result_file=open(s)
sampleLabel=os.path.splitext(s)[0]
events = sum(1 for line in result_file) -1 #quickly determine number of events
result_file.seek(0) #rewind file to beginning
header = result_file.readline()
header = header.strip()
headers = header.split("\t")
pop_offset=len(headers)
popList=[]
for i,header in enumerate(headers):
if header == "pop1":
pop_offset=i
if "pop" in header:
popList.append(header)
markers = headers[0:pop_offset]
result_file.seek(0) #rewind file to beginning
df = pd.read_csv(s, sep='\t')
dataIndex={}
for i,header in enumerate(headers):
dataIndex.update({header:i})
df['pop0']=0
return [sampleLabel,headers,markers,popList,df]
def parseDAFi(s):
df = pd.read_csv(s, sep='\t')
df['pop0']=0
return df
def html_row(file):
return '<img src="{}" style="display:inline;margin:1px" title="{}"/>'.format(export_path+"/"+file+".png",file,file)
display(md("# Automated Analysis Report with Static Composite Dot Plots"))
titlefilename=glob('description.txt')
if titlefilename:
titlefile=open(titlefilename[0])
title=titlefile.readline()
titlefile.close()
display(md("# Dataset: %s"%(title)))
else:
display(md("# Dataset: no description given"))
%%output backend='bokeh'
%%opts Table [width=1200]
metadatafilename=glob('metadata.txt')
if metadatafilename:
metadatafile=open(metadatafilename[0])
metaheader = metadatafile.readline()
metaheader = metaheader.strip()
metaheaders = metaheader.split("\t")
metaDf=pd.read_csv('metadata.txt', sep='\t')
metaTable=hv.Table(metaDf)
display(md("## Metadata"))
display(metaTable)
else:
display(md("## No metadata info given"))
gatedFiles=sorted(glob('Gated/*/flock*.txt'))
gatedDelayed=[[(os.path.split(os.path.dirname(fn))[1]),delayed(parseDAFi)(fn)] for fn in gatedFiles]
sample_labels=[os.path.split(os.path.dirname(fn))[1] for fn in gatedFiles]
dfArray=compute(*gatedDelayed, get=dask.threaded.get)
headers=list(dfArray[0][1])
pop_offset=len(headers)
popList=[]
for i,header in enumerate(headers):
if header == "pop1":
pop_offset=i
if "pop" in header:
popList.append(header)
markers = headers[0:pop_offset]
%%output backend='bokeh'
%%opts Table.gates [width=1200]
%%opts Table.summary [width=1200]
configLabel="pipeline.config"
gates=config_objects(configLabel)
num_gates = len(gates)
summary=config_summary(configLabel, headers)
num_gates = len(summary)
gatesummary = [v for v in summary.values()]
di = {summary.get(element)[0]:str(summary.get(element)[0]).zfill(2)+"_"+summary.get(element)[4] for i,element in enumerate(summary)}
summaryTable=hv.Table(gatesummary,kdims=['Population','Parent','XMarker','YMarker','phenotype','startx', 'endx', 'starty', 'endy'], group='summary', label='Summary')
sortedTable=summaryTable.sort('Population')
sortedTable
%%output backend='bokeh'
axis_popIndexDict = defaultdict(list)
popBounds={}
axises=[]
composite_axis=0
last_xmarker=""
last_ymarker=""
last_parent=0
gatesconfig=[]
for i in range(len(gates)):
pop="pop"+str(i+1)
config=gates.get(pop)
xmarker=str(headers[config[1]-1])
ymarker=str(headers[config[2]-1])
startx=int((float(config[3])/200)*4096)
starty=int((float(config[5])/200)*4096)
endx=int((float(config[4])/200)*4096)
endy=int((float(config[6])/200)*4096)
parent=int(config[7])
ctype=int(config[8])
popBounds.update({pop:[xmarker, ymarker, startx,starty,endx,endy,ctype,"pop"+str(parent)]})
key="axis"+str(composite_axis).zfill(2)
if (xmarker != last_xmarker) or (ymarker != last_ymarker) or (parent != last_parent):
composite_axis=composite_axis+1
key="axis"+str(composite_axis).zfill(2)
axises.append([xmarker,ymarker,key,"pop"+str(parent)])
axis_popIndexDict[key].append(pop)
last_xmarker=xmarker
last_ymarker=ymarker
last_parent=parent
gatesconfig.append([pop,xmarker,ymarker,parent])
num_axises = len(axises)
markerTable=hv.Table(markers,kdims=['Marker'])
axis_popTable=hv.Table(axis_popIndexDict, kdims=['Axis Index'], vdims=['sub populations'])
markerTable+axis_popTable.sort('Axis Index')
hv.notebook_extension('matplotlib')
batchpercent_df = pd.read_csv('Gated/Batch_percentages.txt', sep='\t', index_col=0)
batchpercent_df=batchpercent_df.rename(di).round(2)
with pd.option_context('display.max_columns', None):
display(batchpercent_df)
batchevents_df = pd.read_csv('Gated/Batch_events.txt', sep='\t', index_col=0)
batchevents_df=batchevents_df.rename(di)
with pd.option_context('display.max_columns', None):
display(batchevents_df)
%%output backend="bokeh"
%%opts Table [width=1000]
p_df=pd.DataFrame(batchpercent_df.unstack())
p_df.columns=['Percent']
e_df=pd.DataFrame(batchevents_df.unstack())
e_df.columns=['Events']
c_df=pd.concat([p_df,e_df],axis=1, join='outer').reset_index()
c_df.columns=['Sample','Population','Percent','Events']
c_df=c_df.replace({"Population":di})
%%output backend="bokeh" size=150
%%opts BoxWhisker [xrotation=45]
percentBoxPlot=hv.BoxWhisker(c_df, kdims=['Population'],vdims='Percent').relabel('Population Percent Box Plot')
eventsBoxPlot=hv.BoxWhisker(c_df, kdims=['Population'],vdims='Events').relabel('Population Events Box Plot')
percentBoxPlot+eventsBoxPlot
# centroidFiles=sorted(glob('Gated/*/cent*.txt'))
# reclusterset=set(["pop"+os.path.basename(fn).split(".")[0].split("centroids")[1]for fn in centroidFiles])
# reclustermap={}
# currentcluster="pop0"
# for j, gate in enumerate(gatesconfig):
# parent="pop"+str(gate[3])
# if parent in reclusterset:
# currentcluster=parent
# reclustermap.update({gate[0]: currentcluster})
# centDict=dict(((os.path.split(os.path.dirname(fn))[1]+"_pop"+os.path.basename(fn).split(".")[0].split("centroids")[1], parseDAFi(fn)) for fn in centroidFiles))
# #centDict=compute(*centDelayed, get=dask.threaded.get)
# centArray=[[sample, gate[0], centDict.has_key(sample+"_"+reclustermap.get(gate[0]))] for k, sample in enumerate(sample_labels) for j, gate in enumerate(gatesconfig)]
poplist=natural_sort(popBounds.keys())
hv.opts("RGB [width=600, height=600, bgcolor='#D3D3D3', fontsize={'title':15, 'xlabel':10, 'ylabel':10, 'ticks':10}]")
hv.opts("Points.cent (color='purple' marker='+' size=10)")
size=400
popdfPlots = hv.HoloMap({(sample, j+1): datashade(hv.Points(dfArray[k][1].loc[dfArray[k][1][gate[0]]==0], kdims=[gate[1], gate[2]]), width=size, height=size, x_range=(0,4096), y_range=(0,4096), dynamic=False, link_inputs=False, cmap=cmap_pop[j])
for k, sample in enumerate(sample_labels) for j, gate in enumerate(gatesconfig)}, kdims=['Sample', 'Pop'])
alldfPlots = hv.HoloMap({(sample, j+1): datashade(hv.Points(dfArray[k][1], kdims=[gate[1], gate[2]]), width=size, height=size, x_range=(0,4096), y_range=(0,4096), dynamic=False, link_inputs=False, cmap=cmap_all)
for k, sample in enumerate(sample_labels) for j, gate in enumerate(gatesconfig)}, kdims=['Sample', 'Pop'])
parentdfPlots = hv.HoloMap({(sample, j+1): datashade(hv.Points(dfArray[k][1].loc[(dfArray[k][1]["pop"+str(gate[3])]==0) & (dfArray[k][1][gate[0]]==1)], kdims=[gate[1], gate[2]]), width=size, height=size, x_range=(0,4096), y_range=(0,4096), dynamic=False, link_inputs=False, cmap=cmap_parent)
for k, sample in enumerate(sample_labels) for j, gate in enumerate(gatesconfig)}, kdims=['Sample', 'Pop'])
boundarydfPlots = hv.HoloMap({(sample, j+1): hv.Bounds((popBounds.get(gate[0])[2], popBounds.get(gate[0])[3], popBounds.get(gate[0])[4], popBounds.get(gate[0])[5])).opts(style=dict(line_color=cmap_pop[j][0],color=cmap_pop[j][0]))
for k, sample in enumerate(sample_labels) for j, gate in enumerate(gatesconfig)}, kdims=['Sample', 'Pop'])
# centroiddfPlots = hv.HoloMap({(sample, j+1): hv.Points(centDict.get(sample+"_"+reclustermap.get(gate[0])), kdims=[gate[1], gate[2]], group="cent")
# for k, sample in enumerate(sample_labels) for j, gate in enumerate(gatesconfig)}, kdims=['Sample', 'Pop'])
def outputSampleGates(sample):
filename=export_path+"/"+sample
hv.output(hv.NdLayout(combineddfPlots[sample,:]).cols(1), backend='matplotlib', size=200, fig='png', filename=filename)
#print(filename)
return filename
def outputPopGates(pop):
filename=export_path+"/"+pop
length=len(combineddfPlots[:,poplist.index(pop)])
hv.output(hv.NdLayout(combineddfPlots[:,poplist.index(pop)]).cols(length), backend='matplotlib', size=200, fig='html', filename=filename)
return filename
def displaySampleGates(sample):
filename=export_path+"/"+sample
temp=hv.NdLayout(combineddfPlots[sample,:]).cols(1)
return temp
def displayPopGates(pop):
filename=export_path+"/"+pop
length=len(combineddfPlots[:,poplist.index(pop)])
temp=hv.NdLayout(combineddfPlots[:,poplist.index(pop)]).cols(length)
return temp
#hv.opts("RGB [width=600, height=600, bgcolor='#D3D3D3', fontsize={'title':15, 'xlabel':10, 'ylabel':10, 'ticks':10}]")
#temp=hv.output(hv.NdLayout(combineddfPlots).cols(5), backend='matplotlib', size=200, fig='png', filename=export_path+"/sample_composite")
#display(HTML('<iframe src=PNG/pop2.html width=3000 height=350> </iframe>'),HTML('<iframe src=PNG/pop3.html width=3000 height=350></iframe>'))
#%%time
#testOutput2=[[pop, outputPopGates(pop)] for pop in poplist]
combineddfPlots=alldfPlots*parentdfPlots*popdfPlots*boundarydfPlots
#combineddfPlots=alldfPlots*parentdfPlots*popdfPlots*boundarydfPlots*centroiddfPlots
testOutput=[[sample, outputSampleGates(sample)] for sample in sample_labels]
#html=''.join(html_row(sample) for sample in sample_labels)
#widgets.HTML(value=html, layout=row_layout)
#testOutput2=[[pop, outputPopGates(pop)] for pop in poplist]
#images=[Image(filename = export_path+"/"+pop+".png") for pop in poplist]
#display(*images)
html="".join(html_row(sample) for sample in sample_labels)
html="<div style=\"width:100% !important; overflow-x: auto;white-space: nowrap;\">"+html+"</div>"
display(HTML(html))